-uonoinfc^ 



Mem. S.A.It. Vol. 75, 282 
©SAIt 2008 



Memorie deiia "/ 




"'^OegliSpO'^ 



Improving Large-scale Convection 
Zone-to-Corona Models 

W. p. Abbetti and G. H. Fisher^ 



o 



o 



p 

o 
o 



X 



Space Sciences Laboratory, University of California, Berkeley, CA 94720-7450 e-mail: 
abbett@ssl.berkeley.edu 



Abstract. We introduce two new methods tliat are designed to improve tiie realism and util- 
ity of large, active region-scale 3D MHD models of the solar atmosphere. We apply these 
methods to RADMHD, a code capable of modeling the Sun's upper convection zone, pho- 
tosphere, chromosphere, transition region, and corona within a single computational vol- 
ume. We first present a way to approximate the physics of optically-thick radiative transfer 
without having to take the computationally expensive step of solving the radiative transfer 
equation in detail. We then briefly describe a rudimentary assimilative technique that allows 
a time series of vector magnetograms to be directly incorporated into the MHD system. 
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1. Introduction 

In this paper, we briefly summarize our ef- 
forts to improve our models of quiet Sun 
and active region magnetic fields in computa- 
tional domains that include the upper convec- 
tion zone, photosphere, chromosphere, transi- 
tion region and low corona within a single 
computational do main. Our goa l is similar to 
fliat presented in lAbbetd (l2007h — that is, to 
develop the techniques necessary to efficiently 
simulate the spatially and temporally disparate 
convection zone-to-corona interface over spa- 
tial scales sufficiently large to accommodate at 
least one active region. 

The advantage of this type of single- 
domain modeling is clear. For example, evolv- 
ing a turbulent convection zone and corona 
simultaneously in a physically self-consistent 
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way allows for the quantitative study of impor- 
tant physical processes such as flux emergence, 
submergence and cancellation; the transport of 
magnetic free energy and helicity into the solar 
atmosphere; the generation of magnetic fields 
via a convective dynamo; and the physics of 
coronal heating. 

However, this approach is challenging. The 
computational domain is highly stratified — 
average thermodynamic quantities change by 
many orders of magnitude as the domain 
transitions from a relatively cool, turbulent 
regime below the visible surface, to a hot, 
magnetically-dominated and shock-dominated 
regime high in the model atmosphere. In ad- 
dition, the low atmosphere is where the ra- 
diation field transitions from being optically 
thick to optically thin. The chromosphere itself 
presents an additional challenge, since the ra- 
diation field is often decoupled from the ther- 
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mal pool, particularly in some of the strongest, 
most energetically important transitions. 

There are a number of ways to model the 
energetics of the convection zone-to-corona 
system, ranging from approximate, parameter- 
ized descri2tions_of_tiieJhemiodynamics (see 
e.g.. lFanll2009l: lHood et al."2009'). to highly re- 
ahstic treatments of radiative transfer (see e.g., 
iMartinez-Svkora et al.i 20Q9a b). Since our ob- 
jective is to model the coupled system over 
large spatial scales, our goal is to find the most 
efficient treatment of the energetics possible 
that still provides a physically meaningful rep- 
resentation of the dynamic connection between 
the convection zone and corona. 

In order to describe the thermodynamics 
of the corona, a model should include the 
efifects of electron thermal conduction along 
magnetic field lines and radiative cooling in the 
optically-thin limit. In addition, some physics- 
based or empirically-based source of coronal 
heating must be present if the model corona is 
to remain hot. In the convective interior well 
below the visible surface, radiative cooling can 
be treated in the diffiision limit. The trick is, 
how best to describe the effects of optically- 
thick radiative transfer in the region of the 
model atmosphere that lies between these two 
extremes. 

The most satisfying approach would be to 
couple the LTE transfer equation (or non-LTE 
population and transfer equations) to the MHD 
system to obtain cooling rates and intensities 
that could be compared directly to observa- 
tions. Unfortunately, for large active region or 
global-scale problems, the computational ex- 
pense of these techniques remains prohibitive. 

In lAbbeta (l2007l) we tried the opposite ap- 
proach — ignore the transfer equation alto- 
gether, and develop an artificial, fully param- 
eterized means of approximating surface cool- 
ing (in this case, we employed a modified form 
of Newton cooling). This worked relatively 
well, provided we carefully calibrated the ad- 
justable parameters to match the average sub- 
surface stratification of previous, more realis- 
tic simulations of magneto-convection where 
the LTE transfer equation was solved in detail 
(lBercikll2002) . 



Of course, the principle drawback of this 
approach is that it is ultimately ad hoc and 
unphysical, and requires other, more realistic 
simulations as a basis for calibration in order to 
get meaningful results. We therefore have de- 
veloped an approximation that is based on the 
macroscopic radiative transfer equation, and 
have incorporated this new treatment into our 
3D MHD model, RADMHD. We describe this 
new method in Section 2. 

While it is important to treat the energetics 
of the system in a physically meaningful way, 
it is also important to remember that the utility 
of a given simulation ultimately depends on the 
statement of the problem. For an MHD simu- 
lation, this boils down to one's choice of initial 
states and boundary conditions. It is of great 
benefit, for example, to pose a simple, well- 
defined problem, and set up a numerical exper- 
iment that can shed light on what is believed to 
be the relevant physical processes in an other- 
wise complex system. For example, important 
progress has been made in understanding the 
physics of magnetic flux emergence by study- 
ing how idealized twisted flux ropes emerge 
through h ighly- stratified mod e l atmospheres 
(see e.g., ICheung et all 120071: iFan & GibsonI 
12004 . 

Yet the observed evolution of the photo- 
spheric magnetic field is often far more com- 
plex, particularly in and around CME and flare 
producing active regions. It is very difficult to 
set up a simple magnetic and energetic con- 
figuration that can initialize a simulation that 
will faithfully mimic the observed evolution of 
a real active region. It is desirable to do so, 
however, since we wish to quantitatively un- 
derstand the physical mechanisms of energy 
storage and release, and the transport of mag- 
netic energy and helicity between the convec- 
tive interior and corona. 

To make progress, we could take a cue 
from meteorologists, and investigate a means 
to incorporate observational data directly into 
MHD models. This is not at all straightfor- 
ward for solar models however, since data is 
obtained entirely through remote sensing, and 
not in situ. 

To address this challenge, we have devel- 
oped a simple, rudimentary means of assimilat- 
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ing a time series of vector magnetograms into 
an MHD model of the photosphere-to-corona 
system. We briefly summarize this technique in 
Section 3, and apply it to the specific problem 
of finding a 3D magnetic field that is as force- 
free as possible given a single measurement of 
the vector magnetic field at the photosphere. 

2. An Approximate Treatment of 
Optically Thick Cooling 

What follows is a brief description of our ap- 
proximate treatment of optically-thick radia- 
tive cooling in the portion of the computational 
domain that represents the solar photosphere 
and chromosphere. In practice, this cooling is 
incorporated into the MHD system as a source 
term in the equation that evolves the internal 
energ y per unit volume (Equation 4 of lAbbetj 
l2Q07i) . 

We begin by characterizing the net cooling 
rate for a volume of plasma at some location 
within the solar atmosphere: 



"-hf 



dV (rjy - Kyly) . 



(1) 



Here, v represents frequency, and Q. solid 
angle. The emissivity, opacity, and specific in- 
tensity are frequency dependent, and are de- 
noted rjy, Ky , and ly respectively. Rearranging 
the order of integration, and defining the source 
function Sy as the ratio of the emissivity to 
opacity, we have 



R 



-jdy^yj, 



dQ.{Sy-Iy). 



(2) 



Since the source function is independent of di- 
rection, we recast the integral as 



R = 4n I dVKyiSy - Jy) 

with mean intensity 

Jy = — fd^Iy. 

An J 

The formal solution for the specific intensity in 
the plane-parallel approximation is 
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Fig. 1. Shown is a comparison of the average 
temperature stratification between the reali stic ra- 
diative magneto-convection simulations of BercikI 
( l2002h (crosses), and the RADMHD model convec- 
tion zone using the new treatment of optically-thick 
transfer (diamonds). 



where fi is the usual cosine angle. Then the 
mean intensity can be expressed as 
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The integral over /i can now be evaluated and 
the mean intensity can be cast as 
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1 r" 

2 Jo 



dT'Sy(T')Ei(\Ty-T'\), 



(7) 



where Ei denotes the first exponential integral. 
So far, th is is simply tex tbook radiative trans- 
fer (e.g., lMihalaslll978h . No approximations 
have yet to be made, other than an assumption 
of a locally plane-parallel atmosphere. Now 
we'll make our first approximation. Note that 
E\(\Ty - t'\) is singular when r' = Ty, and that 
the singularity is integrable. Since Ei is peaked 
around t,,, contributions from Syir') will be 
centered around Sy(Ty). Thus, we approximate 
the mean intensity by 



Jy 



1 



Jo 



Sy(Ty) dT'Eli\Ty-T'\). 



(8) 



This integral can then be evaluated, giving 
a simple expression for the mean intensity. 



(5) Jy~SyiTy)\l- 



EliTy) 



(9) 
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Fig. 2. The temperature along a horizontal slice through a RADMHD model photosphere during the relax- 
ation process. This simulation uses the ad hoc Newton cooling described in Abbett (2007.) to approximate 
optically-thick surface cooling. The slice spans 75 x 37.5 Mm^. 



where £2 refers to the second exponential inte- 
gral. Note that this can be rewritten in the fol- 
lowing way: 



7v _ E2{Ty) 



(10) 



We now return to our expression for the net 
cooling rate and recast it in a slightly different 
form. 



R = 47: 



I dvK^ 



SJl 



Jy 



(11) 



Substituting equation [TO]into the integrand, we 
have 



R ^2n \ dv K^ 



SyEziTy). 



(12) 



If we further assume LTE, the source function 
is simply the Planck function and we have 



R ^ 2n I dvK^ 



By(T)E2iTy). 



(13) 



Now for the real swindle! Let's integrate 
over frequency, and replace the frequency- 
dependent opacity by its Planck weighted 



mean value: That is, replace KyBy with 
K{o-/7r)T'^ where k represents a Planck- 
weighted mean opacity, and cr the Stefan- 
Boltzmann constant. Then including the expo- 
nential function in the average over frequency, 
we find that 



R^2CKO-T'^E2{aT). 



(14) 



Here, C represents the normalization constant 
for the integration. The arbitrary constant a ap- 
pears in the exponential integral since the mean 
opacity used in the calculation of the optical 
depth scale could differ in general from the 
mean opacity that appears by itself in the in- 
tegrand. 

To determine the normalization constant C, 
we integrate our cooling function from zero to 
infinity over an isothermal slab to obtain the 
total radiative flux. The resulting expression 
must be equal to the known result F,o, = crT'^. 
This allows us to determine the normalization 
constant C = a. To evaluate a, we compare 
the detailed cooling rate depth distribution us- 
ing this formu lation with the cooling rate in the 
iBercikI ( |2002 |) LTE model of the solar atmo- 
sphere, and conclude the best-fit value is a = 1 
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Fig. 3. The temperature along a horizontal slice through a RADMHD model photosphere during the relax- 
ation process. In this case, we use our new treatment to approximate the optically-thick surface cooling. We 
note that the model convection zone has yet to fully relax. 



(See Figure[TJ. Thus, our approximate cooling 
function takes the form: 



Ra2Ko-T^E2(T). 



(15) 



The advantage of this treatment lies in its 
simplicity. The above approximation for sur- 
face cooling, while non-linear, is trivial to cal- 
culate for each mesh element. It is certainly 
more phy sical than the ad hoc treatment em- 
ployed in lAbbetd (l2007h . since it is based on 
the radiative transfer equation, and incorpo- 
rates an optical depth scale into the model. 
Further, it has no adjustable parameters. The 
only calibration now required is a choice of 
optical depth ranges over which to apply the 
different approximations. Currently, we use the 
radiative diffusion approximation for optical 
depths greater than 10, an expression for radia- 
tive cooling in the optically-thin limit for op- 
tical depths less than 0.1, and the above treat- 
ment in the intervening layers. 

Figures |2] and [3] provide a qualitative com- 
parison between a model convection zo ne gen- 
erated using the ad hoc approach of lAbbetJ 
(l2007h to estimate the effects of optically-thick 
radiative surface cooling, and one that utilizes 
the approximation described above. The im- 
ages represent the gas temperature along a 



horizontal slice through a RADMHD model 
photosphere during the relaxation process. 
Both simulations used the same initial convec- 
tive state and boundary conditions (periodic 
in the horizontal directions and closed verti- 
cally); the only difference is the treatment of 
the optically-thick transfer Distinct differences 
rapidly develop — the convective cells become 
more irregularly shaped while the size distri- 
bution of cells begins to more closely mimic 
that of the more realistic simulations of iBercikI 
(I2002h . 

However, there are irregularities in the cur- 
rent data set. For example, there are regions 
within the intergranular lanes that are hotter 
than expected. This may be an artifact result- 
ing from our empirically-ba sed coronal heating 
function (see lAbbetJl2007l for details) extend- 
ing unphysically deep into the atmosphere (i.e., 
its optical depth cutoff is too high), or it may 
simply be a transient effect that will subside 
as the simulation progresses. This is a work 
in progress, and we continue to test and val- 
idate the new treatment against our previous 
quiet Sun simulations, and against more real- 
istic magnetoconvection simulations that treat 
the LTE transfer equation in detail. 
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3. Rudimentary Data Assimilation 

We now turn our attention to the problem of 
incorporating a time series of vector magnetic 
field measurements into our MHD model of 
the solar atmosphere. The essential problem 
is that even the most carefully pre-processed 
sequences of vector magnetograms cannot be 
expected to exactly satisfy Faraday's law, and 
thus are physically inconsistent from the point 
of view of the numerical model. 

This is not particularly surprising, since 
it is a non-trivial task to properly transform 
polarization measurements into a vector mag- 
netic field. The datasets naturally suffer from 
the effects of uncertainty due to noise, see- 
ing, or saturation; the inversion process itself is 
model dependent; and the well-known 180 de- 
gree ambiguity in the transverse field must be 
resolved in the context of a timeseries of mag- 
netograms rather than in a single magnetogram 
in isolation. 

While the data itself presents challenges, it 
is important to remember that the model suffers 
from fundamental deficiencies of its own. The 
single-fluid MHD system may not capture the 
essential physics of the solar atmosphere, par- 
ticularly in the low atmosphere where effects 
such as non-LTE transfer, ion-neutral diffu- 
sion, magnetic reconnection, and non-thermal 
partical acceleration may play fundamental 
roles in the dynamic evolution of a given re- 
gion. 

Then how shall we proceed? One approach 
is to take the data at face value and incorpo- 
rate the measurements directly into an MHD 
model in the form of time-de pendent, charac - 
teristic boundary conditions (IWu et al.ll2006h . 
Here, one must be mindful to not over-specify 
the MHD system. This method is restrictive 
in that only certain components of the electric 
field or flow inferred from the data can be used 
to drive the simulation. In addition, one must 
make assumptions about the thermodynamics 
of the system in order to drive the model atmo- 
sphere in a physical way. 

Here, we take another approach. To avoid 
the mathematical constraints inherent to MHD 
boundary conditions, we push our lower 
boundary slightly deeper into the photosphere 



and instead incorporate the data into the model 
via additional forces acting on active zones of 
the calculation where the entire MHD system 
is being self-consistently evolved. 

To do this, we must first obtain an in- 
ductive flow field from a given time series 
of magnetograms that is both consistent with 
the observed evolution of the vector field and 
Faraday's law. This is a non-trivial task, as 
the problem is inherently under-determined, 
and the cadence and quality of the magne- 
tograms may vary. There is a growing num- 
ber of inversion techniques that address this 

I 1 r — 1 I 

problem (Fisher et al. 2009; Ravindra et al. 
2008; Schuck 2008; Georgoulis & LaBonte 



2006; Welschetal. 2004t iLongcopd 12004 : 
iKusano et al.i i2002) ; each is capable of pro- 
viding an inductive flowfield consistent with 
the observations and suitable for incorporation 
into an MHD model. 

Next, we must generate an initial near- 
equilibrium or steady state atmosphere from 
the first magnetogram of the timeseries, and 
choose an appropriate set of exterior bound- 
ary conditions. The challenge here is to mini- 
mize perpendicular currents within the compu- 
tational volume, and to provide coronal bound- 
ary conditions that minimize forces resulting 
from magnetic tension. This way, the forces 
introduced into the model photosphere are the 
principle drivers of the system. 

As a starting point, we generate a non- 
constant-a force-free extrapolation using a 
variation of the opti mization technique of 
IWheatland et all (12000). Given the photo- 
spheric magnetogram and a choice of external 
boundary conditions, this procedure minimizes 
the functional 



/-( 



|(VxB)xBp 



H-IV-BI 



(16) 



and generates an initial magnetic configura- 
tion. 

Unfortunately, the reality is that the pho- 
tosphere is often far from force-free, making 
the mathematical problem of generating per- 
fectly force-free equilibia ill-posed. While the 
optimization techniqu e performs well relat ive 
to other methods (see ISchrijver et al.ll2()()6h . it 
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Fig. 4. The three components of the magnetic field (from left to right B-, S „ and By) from the early stages 
of a 3D RADMHD simulation of NOAA AR 8210. The top row shows a horizontal slice through the model 
photosphere, and the bottom row shows the magnetic field through a horizontal slice corresponding to the 
model chromosphere. The top row mirrors the vector magnetic field as measured by the IVM vector magne- 
tograph at Mees Solar Observatory on Haleakala, and the bottom row represents a synthetic chromospheric 
magnetogram as predicted by the MHD model. 



still cannot be expected to fully converge to 
an equilibrium state without altering the trans- 
verse magnetic field at the photospheric bound- 
ary. 

We therefore use the optimization method 
to generate an initial starting point for an MHD 
relaxation. The above functional need not be 
vanishingly small in every mesh element, since 
the MHD code will difiiise away any signifi- 
cant divergence error, and clean up any noisy, 
unphysical currents near the lower boundary 
(see Figure nil. In practice, this is done by ar- 
tificially damping fast-moving waves and al- 
lowing the system to slowly evolve to a near- 
equilibrium state. Of course, the resulting at- 
mosphere is not expected to be force-free near 
the photosphere. The currents in the system 
are, however, more physical since they were 
evolved via the MHD system of conservation 



equations rather than by attempting to mini- 
mize the functional of equation[T6l 

Our purpose here is not to find a perfectly 
stable equilibrium solution. In fact, such a state 
may not exist, given the vector magnetogram 
and choice of boundary conditions. We are 
simply striving for an initial atmosphere that 
is not so vastly out of force balance that mo- 
tions at the model photosphere are immedi- 
ately overwhelmed by other less relevant pro- 
cesses. Once this is achieved, we drive the at- 
mosphere in the following way. 

First, we define the physical contribution 
to the force as that described by the MHD 
momentum conservation equation (see lAbbetj 
|2007| for details), 



Fs -V- 



puu + (;, + -I---n 



+pg (17) 
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We then define the forces impHed by the data, 

dpUiny 



'-data — 



dt 



(18) 



Here, uinv refers to the inductive flow field ob- 
tained through one of the many velocity inver- 
sion techniques (see lWelsch et alfcOOTI) . 

Then in a thin volume corresponding to the 
model's photosphere, we recast the momentum 
equation in the following form: 

^1 , = ^i^dam)^ + (1 - ^) (F)^ + (F)|| (19) 
Ot 'phot 

where the parallel and perpendicular subscripts 
denote the forces parallel or perpendicular to 
the direction of the magnetic field. Here, < 
^ < 1 represents a "confidence matrix" de- 
fined at each mesh element within the pho- 
tospheric volume. It is easy to see that when 
^ = 1, the forces perpendicular to the magnetic 
field within the model photosphere are deter- 
mined entirely by the data, and when ^ = 0, the 
photospheric layer evolves as the MHD system 
normally would in the absence of any observa- 
tional forcing. 

Since flows parallel to the field do not 
affect magnetic evolution, we allow them to 
evolve in an unconstrained fashion. All other 
independent variables evolve as prescribed by 
the MHD system of equations including the 
magnetic field. Recall, uinv is designed to be 
consistent with of the observed evolution of 
one or more components of the photospheric 
magnetic field (depending on the inversion 
method used) and, in principle, should drive 
the photospheric field in a manner consistent 
with the timeseries of magnetograms. 

4. Concluding Remarks 

We have developed a rudimentary means of 
assimilating a time series of vector magne- 
tograms into the interior volume of an MHD 
model in a manner that is stable, and does not 
over-specify the problem. We are currently us- 
ing this assimilative technique to incorporate 
a timeseries of vector magnetograms into a 
120 Mm^ RADMHD model atmosphere that 
contains a model photosphere, chromosphere. 



transition region and corona. The IVM data we 
are using is a four hour timeseries from NOAA 
AR 8210 — a wefl-studied flare and CME- 
producing active region. The simulations are in 
their preliminary stages, and we hope to report 
on this work in the near future. 

In addition, we have presented a com- 
putationally efficient method of approximat- 
ing optically-thick radiative cooling in our 
RADMHD quiet Sun models. The treatment 
improves upon the method of I AbbettI (l2007h . 
while still retaining the efficiency necessary to 
allow for large, active region-scale, convection 
zone-to-corona computational domains. Our 
simulations are progressing, and we are cur- 
rently evaluating the efficacy and reliability of 
the new method. We are optimistic that each 
of these methods will improve the realism and 
utility of our current suite of numerical models. 
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